program principal
interface

subroutine trapecioCompuestoFichero(integVG1,integVG2,vector,VG1,VG2,n)
	real(kind=8),dimension(:),allocatable::vector,VG1,VG2	
	real(kind=8)::integVG1,integVG2
	integer::i,n
end subroutine

subroutine simpsonCompuestoFichero(integVG1,integVG2,vector,VG1,VG2,n)
	real(kind=8),dimension(:),allocatable::vector,VG1,VG2
	real(kind=8)::integVG1,integVG2
	integer::i,n
end subroutine

subroutine trapecioCompuesto(inter,integ,a,b)
	real(kind=8)::a,b,valora,valorb,integ
	integer::inter
end subroutine

subroutine simpsonCompuesto(inter,integ,a,b)
	real(kind=8)::a,b,valora,valorb,valorc,integ
	integer::inter
end subroutine


function trapecio(a,b,valora,valorb)result(res3)
	real(kind=8)::a,b,valora,valorb,res3
end function

function simpson(a,b,valora,valorc,valorb)result(res4)
	real(kind=8)::a,b,valora,valorb,valorc,res4
end function

subroutine datos(a,b,maxit,tol,formula,tec,inter)
	real(kind=8)::tol,a,b
	integer::maxit,formula,tec,inter
end subroutine

function funcion(x,integVG1,integVG2)result(res1)
	real(Kind=8),intent(in)::x,integVG1,integVG2
	real(kind=8)::res1
end function


subroutine dicotomia(a,b,it,tol,integVG1,integVG2)
	real(kind=8)::tol,a,b,x,ant,integVG1,integVG2
	integer::it,i
end subroutine

function LeerDatos(v,valores,ea)result(n)
	real(kind=8),dimension(:),allocatable::v,valores
	integer::n
	real(kind=8)::ea
end function

subroutine Lagrange(x,w,dimen,matriz)
	real(kind=8),dimension(:),allocatable::x,w
	real(kind=8),dimension(:,:),allocatable::matriz
	integer::dimen
end subroutine

function evaluar(matriz,x,pto,dimen)result(res)
	real(kind=8),dimension(:,:),allocatable::matriz
	real(kind=8),dimension(:),allocatable::x
	real(kind=8)::pto,sol,res
	integer::i,j,dimen
end function

subroutine calcVectAplTrap(vector,imagenes,VG1,VG2,ea,n)
	real(kind=8),dimension(:),allocatable::vector,imagenes,VG1,VG2
	real(kind=8)::ea,pi,radian
	integer::i
end subroutine

subroutine calcVectAplSimp(vector,vinter,imagenes,imaginter,VG1,VG2,ea,n)
	real(kind=8),dimension(:),allocatable::vector,vinter,imagenes,imaginter,VG1,VG2
	real(kind=8)::ea,pi,radian
	integer::i
end subroutine

end interface

real(kind=8),dimension(:),allocatable::v,v2,evaluar2,VG1,VG2,valores
real(kind=8),dimension(:,:),allocatable::matriz
real(kind=8)::integ,integ1,error,a,b,ea,integVG1,integVG2,valora,valorb,valorc,d,evpto
integer::it,formula,dimen,i,inter,tec,n
	
	call datos(a,b,it,error,formula,tec,inter)
	dimen=LeerDatos(v,valores,ea)
	integ=0._8
	integVG1=0._8
	integVG2=0._8
	allocate(v2(0:dimen-1))
	allocate(evaluar2(0:dimen-1))
	if(formula==1)then!simpson
		if(tec==1)then
			call simpsonCompuesto(inter,integ,a,b)
		else
			call Lagrange(v,valores,dimen,matriz)
			do i=0,dimen-1,1	
				v2(i)=(v(i)+v(i+1))*0.5_8
				evpto= evaluar(matriz,v,v2(i),dimen)
				evaluar2(i)=evpto
			end do
			call calcVectAplSimp(v,v2,valores,evaluar2,VG1,VG2,ea,dimen)
			call SimpsonCompuestoFichero(integVG1,integVG2,v,VG1,VG2,(dimen*2))
			call dicotomia(v(0),v(10),it,error,integVG1,integVG2)
		end if
	else!trapecio
		
		if (tec==1)then
			call trapecioCompuesto(inter,integ,a,b)
		else
			call calcVectAplTrap(v,valores,VG1,VG2,ea,dimen)
			call trapecioCompuestoFichero(integVG1,integVG2,v,VG1,VG2,dimen)
			call dicotomia(v(0),v(10),it,error,integVG1,integVG2)
		end if
	end if

end program
	
	


 recursive subroutine datos(a,b,maxit,tol,formula,tec,inter)
	real(kind=8)::tol,a,b
	integer::maxit,formula,tec,inter
	print*,'¿datos por teclado(1) o fichero(2)?'
	read(5,*)tec
	if(tec==1)then
		print*,'¿cuantos intervalos?'
		read(5,*)inter
		print*,'¿cual es el intervalo de partida?'
		read(5,*)a
		read(5,*)b
		if(funcion(a)==0._8) then 
			print*,'la raiz es: ',a
			stop('')
		else	if (funcion(b)==0._8) then 
			print*,'la raiz es: ',b
			stop
		end if
	end if
	print*,'¿numero maximo de iteraciones?'
	read(5,*)maxit
	print*,'¿tolerancia de error?'
	read(5,*)tol
	print*,'¿que formula de cuadratura? simpson(1) o trapecio(2)?'
	read(5,*)formula

	
end subroutine



subroutine trapecioCompuestoFichero(integVG1,integVG2,vector,VG1,VG2,n)
	real(kind=8),dimension(:),allocatable::vector,VG1,VG2	
	real(kind=8)::integVG1,integVG2
	integer::i,n
	do i=0,n-2,1
		integVG1=integVG1 +trapecio(vector(i),vector(i+1),VG1(i),VG1(i+1))
		integVG2=integVG2 +trapecio(vector(i),vector(i+1),VG2(i),VG2(i+1))
	end do
end subroutine

subroutine simpsonCompuestoFichero(integVG1,integVG2,vector,VG1,VG2,n)
	real(kind=8),dimension(:),allocatable::vector,VG1,VG2
	real(kind=8)::integVG1,integVG2
	integer::i,n
	do i=0,n-4,2
		integVG1=integVG1 +simpson(vector(i/2),vector(i/2+1),VG1(i),VG1(i+1),VG1(i+2))
		integVG2=integVG2 +simpson(vector(i/2),vector(i/2+1),VG2(i),VG2(i+1),VG2(i+2))
	end do
	
end subroutine


subroutine trapecioCompuesto(inter,integ,a,b)
	real(kind=8)::a,b,valora,valorb,integ,d
	integer::inter,i
	d=(b-a)/inter
	do i=1,inter,1
		valora=funcion(a)
		b=a+d
		valorb=funcion(b)
		integ=integ +trapecio(a,b,valora,valorb)
		a=b
	end do
end subroutine

subroutine simpsonCompuesto(inter,integ,a,b)
	real(kind=8)::a,b,valora,valorb,valorc,integ,d
	integer::inter,i
	d=(b-a)/inter
	do i=1,inter,1
		valora=funcion(a)
		b=a+d
		valorb=funcion(b)
		valorc=funcion((a+b)*0.5_8)
		integ=integ +simpson(a,b,valora,valorc,valorb)
		a=b
	end do
end subroutine


function trapecio(a,b,valora,valorb)result(res3)
	real(kind=8)::a,b,valora,valorb,res3
	res3=(b-a)*0.5_8*(valora+valorb)
end function

function simpson(a,b,valora,valorc,valorb)result(res4)
	real(kind=8)::a,b,valora,valorb,valorc,res4

	res4=(b-a)/6._8*(valora+4._8*valorc+valorb)
end function

function funcion(x,integVG1,integVG2)result(res1)
	real(Kind=8),intent(in)::x,integVG1,integVG2
	real(kind=8)::res1
	res1=integVG1-x*integVG2
end function

subroutine dicotomia(a,b,it,tol,integVG1,integVG2)
	real(kind=8)::tol,a,b,x,ant,integVG1,integVG2
	integer::it,i
	interface
		function funcion(x,integVG1,integVG2)result(res1)
			real(Kind=8),intent(in)::x,integVG1,integVG2
			real(kind=8)::res1
		end function
	end interface
	if(funcion(a,integVG1,integVG2)*funcion(b,integVG1,integVG2)>0._8) then
		print*,'				el intervalo no es valido'
	else
		open(unit=1, file='integracion.dat',status='unknown')
		x=a
		do i=1,it,1
			ant=x
			x=(a+b)*0.5_8
		
			if(funcion(a,integVG1,integVG2)*funcion(x,integVG1,integVG2)<0) then
				b=x
			else
				a=x
			end if
			error=abs((x-ant)/x)
	
			write(1,*)'iteracion: ',i
			write(1,*)'Aproximacion'
			write(1,'(f17.14)'),x
			write(1,*)'error relativo: '
			write(1,'(e21.14)')error
			write(1,*)
			if(funcion(x,integVG1,integVG2)==0._8 .or. error<tol) then
				exit
			end if
		end do
		close(1)
	end if
end subroutine

function LeerDatos(v,valores,ea)result(n)
	real(kind=8),dimension(:),allocatable::v,valores
	integer::n
	real(kind=8)::ea
	open(unit=3, file='fichero.dat',status='unknown')
	read(3,*)ea
	read(3,*)n
	allocate(v(0:n))
	allocate(valores(0:n))
	read(3,*)v
	read(3,*)valores
end function

subroutine Lagrange(x,w,dimen,matriz)
	real(kind=8),dimension(:),allocatable::x,w
	real(kind=8),dimension(:,:),allocatable::matriz
	integer::dimen,fila,columna
	allocate(matriz(0:dimen,0:dimen))
	matriz=0
	matriz(0:dimen,0)=w
	do columna=1,dimen,1
		do fila=0,dimen-columna,1
			matriz(fila,columna)=(matriz(fila+1,columna-1)-matriz(fila,columna-1))/(x(fila+columna)-x(fila))
		end do
	end do
	do columna=0,dimen,1
		print*,matriz(0,columna)
	end do

	open(unit=2, file='Lagrange.dat',status='unknown')
	do i=0,dimen,1
		write(2,'(f19.14)')matriz(i,0:dimen)
	end do
	close(2)
end subroutine

function evaluar(matriz,x,pto,dimen)result(res)
	real(kind=8),dimension(:,:),allocatable::matriz
	real(kind=8),dimension(:),allocatable::x
	real(kind=8)::pto,sol,m,res
	integer::i,j,dimen
	sol=matriz(0,0)
	m=1._8
	do i=1,dimen,1
		m=(pto-x(i-1))*m
		sol=sol+matriz(0,i)*m
	end do
	res=sol
end function

subroutine calcVectAplTrap(vector,imagenes,VG1,VG2,ea,n)
	real(kind=8),dimension(:),allocatable::vector,imagenes,VG1,VG2
	real(kind=8)::ea,pi,radan
	integer::i,n
	pi=acos(-1.0)
	radian=pi/180._8
	allocate(VG1(0:n))
	allocate(VG2(0:n))
	do i=0,n,1
		VG1(i)=500._8*abs(ea)*imagenes(i)
		VG2(i)=30._8*cos(vector(i)*radian)
	end do
end subroutine


subroutine calcVectAplSimp(vector,vinter,imagenes,imaginter,VG1,VG2,ea,n)
	real(kind=8),dimension(:),allocatable::vector,vinter,imagenes,imaginter,VG1,VG2
	real(kind=8)::ea,pi,radian
	integer::i
	pi=acos(-1.0)
	radian=pi/180._8
	allocate(VG1(0:2*n))
	allocate(VG2(0:2*n))

	do i=0,2*n-2,2
		VG1(i)=500._8*abs(ea)*imagenes(i/2)
		VG1(i+1)=500._8*abs(ea)*imaginter(i/2)
		VG2(i)=30._8*cos(vector(i/2)*radian)
		VG2(i+1)=30._8*cos(vinter(i/2)*radian)
	end do
		VG1(2*n)=500._8*abs(ea)*imagenes(n)
		VG2(2*n)=30._8*cos(vector(n)*radian)
end subroutine







